library(tidyverse)
library(ggpubr)
library(here)
library(patchwork)
temperature_range_data <- read_csv(here("1-data/TempChoice_abundance_dryad.csv"))
Temperature_ranges_plot <- ggplot(data = temperature_range_data, aes(x = day, y = mean_density, group = population, color = as.factor(temperature))) +
geom_point() +
geom_line() +
scale_y_log10() +
xlab("Day") +
ylab("Population density (cells/ml)") +
scale_color_manual(name = "Temperature (°C)",
values = c('#053061', '#2166ac', '#4393c3', '#92c5de',
'#d1e5f0','#fddbc7','#f4a582','#d6604d','#b2182b', '#67001f')) +
theme(text = element_text(size = 15),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
legend.position = "bottom",
plot.margin = unit(c(1,1,1.5,1.2),"cm"))
Temperature_ranges_plot
# load the abundance data
abundance_data <- read_csv(here("1-data/TempAdapt_abundances_dryad.csv"))
# load the morphology data
morphology_data <- read_csv(here("1-data/TempAdapt_morphology_dryad.csv"))
# total number of cells
unique(morphology_data$id) %>% length()
## [1] 94344
# mean number of cells plus standard error
cells_summary <- morphology_data %>% group_by(population, date) %>% summarise(measured_cells = n())
cells_summary %>% ungroup() %>% summarise(mean_measured_cells = mean(measured_cells), se_measured_cells = sd(measured_cells)/sqrt(n()))
## # A tibble: 1 x 2
## mean_measured_cells se_measured_cells
## <dbl> <dbl>
## 1 513. 26.9
# calculate mean, variance and sd of area and AR for each population for each day
morphology_data_ts <- morphology_data %>% group_by(population, lineage, batch, date, day, temperature) %>%
summarise(mean_area_day = mean(mean_area), sd_area_day = sd(mean_area),
se_area_day = sd(mean_area)/sqrt(n()), var_area_day = var(mean_area),
mean_AR_day = mean(mean_ar), sd_AR_day = sd(mean_ar),
se_AR_day = sd(mean_ar)/sqrt(n()))
# calculate coefficient of variation of cell area
morphology_data_ts <- morphology_data_ts %>% mutate(cv_area_day = sd_area_day/mean_area_day)
# merge both datasets
morphology_abundance_data <- left_join(abundance_data, morphology_data_ts)
# add name column for plotting (lineage + temperature)
morphology_abundance_data <- morphology_abundance_data %>% mutate(name = paste0(lineage, "_", temperature))
# number of generations per batch ####
generations <- morphology_abundance_data %>% group_by(population, lineage, batch,replicate, temperature) %>%
summarise(maxi_density = max(mean_density), min_density = min(mean_density),
number_generations = (log(max(mean_density)) - log(min(mean_density)))/log(2))
generations_summary <- generations %>% group_by(batch, temperature) %>% summarise(mean_generations_batch = mean(number_generations))
# create a tibble with the number of generations per batch at each temperature
generations_summary <- generations_summary %>% group_by(batch) %>% pivot_wider(names_from = temperature, values_from = mean_generations_batch, names_prefix = "generations_")
# round the number of generations
generations_summary <- generations_summary %>% mutate(text_plot = case_when(
batch == 1 ~ paste0(round(generations_20, digits = 0), " gen."),
batch == 2 | batch == 3 ~ paste0(round(generations_38, digits = 0), " gen."),
batch > 3 ~ paste0(round(generations_20, digits = 0), " gen. 20°C \n", round(generations_38, digits = 0), " gen. 38°C")))
# define size of small and large fonts for the plot
small_font <- 10
large_font <- 13
plot_abundances <- ggplot(data = morphology_abundance_data) +
geom_point(aes(x = batch_day, y = mean_density, group = population, colour = name), size = 1) +
geom_line(aes(x = batch_day, y = mean_density, group = population, colour = name)) +
geom_label(data = generations_summary, aes(x = 9, y = 500, label = text_plot), size = 2.1) +
scale_y_log10() +
ylab("Population density\n(cells/ml)") +
xlab("") +
scale_color_manual(values = c("#d2e0ee", "#efd0d4", "#90b2d5", "#d88b95", "#4d84bc","#c14655", "#2166ac","#b2182b", "red", "blue"),
guide = F, name = "Population and\n temperature",
labels = as_labeller(c("1_20" = "Pop. 1 20°C", "1_38" = "Pop. 1 38°C",
"2_20" = "Pop. 2 20°C", "2_38" = "Pop. 2 38°C",
"3_20" = "Pop. 3 20°C", "3_38" = "Pop. 3 38°C",
"4_20" = "Pop. 4 20°C", "4_38" = "Pop. 4 38°C"))) +
scale_x_continuous(breaks = c(5, 10)) +
theme(text = element_text(size = large_font),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
legend.position = "bottom",
plot.margin = unit(c(0, 1.5, 0 , 1.5 ),"cm")) +
facet_wrap(~ batch, ncol = 5, labeller = as_labeller(c("1" = "Batch 1", "2" = "Batch 2",
"3" = "Batch 3", "4" = "Batch 4", "5" = "Batch 5")))
# minimum and maximum cell area in the control
area_min <- morphology_abundance_data %>% ungroup() %>% filter(batch == 1 & day > 0) %>%
select(mean_area_day) %>% min()
area_max <- morphology_abundance_data %>% ungroup() %>% filter(batch == 1 & day > 0) %>%
select(mean_area_day) %>% max()
# plot cell area
plot_cell_area <- ggplot(data = morphology_abundance_data, aes(x = batch_day, y = mean_area_day,
group = population, col = name)) +
geom_point(size = 1) +
geom_hline(aes(yintercept = area_min), linetype = "dashed") +
geom_hline(aes(yintercept = area_max), linetype = "dashed") +
geom_errorbar(aes(ymin = mean_area_day - se_area_day, ymax = mean_area_day + se_area_day, width = 0.9), size = 0.2) +
geom_line() +
xlab("") +
ylab("Cell size (um2)") +
scale_color_manual(values = c("#d2e0ee", "#efd0d4", "#90b2d5", "#d88b95", "#4d84bc","#c14655", "#2166ac","#b2182b"),
name = "Population and\n temperature",
labels = as_labeller(c("1_20" = "Pop. 1 20°C", "1_38" = "Pop. 1 38°C",
"2_20" = "Pop. 2 20°C", "2_38" = "Pop. 2 38°C",
"3_20" = "Pop. 3 20°C", "3_38" = "Pop. 3 38°C",
"4_20" = "Pop. 4 20°C", "4_38" = "Pop. 4 38°C"))) +
scale_x_continuous(breaks = c(5, 10)) +
theme(text = element_text(size = large_font), legend.text = element_text(size = small_font),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
legend.position = "bottom", strip.text = element_blank(),
plot.margin = unit(c(0, 1.5, 0 ,1.5 ),"cm")) +
facet_wrap(~ batch, ncol = 5, labeller = as_labeller(c("1" = "Batch 1", "2" = "Batch 2",
"3" = "Batch 3", "4" = "Batch 4", "5" = "Batch 5")))
# minimum and maximum cv of cell area in the control
cv_area_min <- morphology_abundance_data %>% ungroup() %>% filter(batch == 1 & day > 0) %>%
select(cv_area_day) %>% min()
cv_area_max <- morphology_abundance_data %>% ungroup() %>% filter(batch == 1 & day > 0) %>%
select(cv_area_day) %>% max()
# plot coefficient of variation of cell area
# multiply cv per 100 to obtain percentage
plot_cv_cell_area <- ggplot(data = morphology_abundance_data, aes(x = batch_day, y = cv_area_day*100,
group = population, col = name)) +
geom_point(size = 1) +
geom_hline(aes(yintercept = cv_area_min*100), linetype = "dashed") +
geom_hline(aes(yintercept = cv_area_max*100), linetype = "dashed") +
geom_line() +
xlab("") +
ylab("Coefficient of variation\n of cell size (%)") +
scale_color_manual(values = c("#d2e0ee", "#efd0d4", "#90b2d5", "#d88b95", "#4d84bc","#c14655", "#2166ac","#b2182b"),
name = "Population and\n temperature",
labels = as_labeller(c("1_20" = "Pop. 1 20°C", "1_38" = "Pop. 1 38°C",
"2_20" = "Pop. 2 20°C", "2_38" = "Pop. 2 38°C",
"3_20" = "Pop. 3 20°C", "3_38" = "Pop. 3 38°C",
"4_20" = "Pop. 4 20°C", "4_38" = "Pop. 4 38°C"))) +
scale_x_continuous(breaks = c(5, 10)) +
scale_y_continuous(breaks = c(15, 30, 45, 60)) +
theme(text = element_text(size = large_font), legend.text = element_text(size = small_font),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
legend.position = "bottom", strip.text = element_blank(),
plot.margin = unit(c(0, 1.5, 0, 1.5 ),"cm")) +
facet_wrap(~ batch, ncol = 5, labeller = as_labeller(c("1" = "Batch 1", "2" = "Batch 2",
"3" = "Batch 3", "4" = "Batch 4", "5" = "Batch 5")))
# minimum and maximum variance of cell area in the control
var_area_min <- morphology_abundance_data %>% ungroup() %>% filter(batch == 1 & day > 0) %>%
select(var_area_day) %>% min()
var_area_max <- morphology_abundance_data %>% ungroup() %>% filter(batch == 1 & day > 0) %>%
select(var_area_day) %>% max()
# plot variance of cell area
plot_var_cell_area <- ggplot(data = morphology_abundance_data, aes(x = batch_day, y = var_area_day, group = population, col = name)) +
geom_point(size = 1) +
geom_hline(aes(yintercept = var_area_min), linetype = "dashed") +
geom_hline(aes(yintercept = var_area_max), linetype = "dashed") +
geom_line() +
xlab("Day") +
ylab("Variance of\n cell size") +
scale_color_manual(values = c("#d2e0ee", "#efd0d4", "#90b2d5", "#d88b95", "#4d84bc","#c14655", "#2166ac","#b2182b"),
name = "Population and\n temperature",
labels = as_labeller(c("1_20" = "Pop. 1 20°C", "1_38" = "Pop. 1 38°C",
"2_20" = "Pop. 2 20°C", "2_38" = "Pop. 2 38°C",
"3_20" = "Pop. 3 20°C", "3_38" = "Pop. 3 38°C",
"4_20" = "Pop. 4 20°C", "4_38" = "Pop. 4 38°C"))) +
scale_x_continuous(breaks = c(5, 10)) +
theme(text = element_text(size = large_font), legend.text = element_text(size = small_font),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
legend.position = "bottom", strip.text = element_blank(),
plot.margin = unit(c(0, 1.5, 0 ,1.5 ),"cm")) +
facet_wrap(~ batch, ncol = 5, labeller = as_labeller(c("1" = "Batch 1", "2" = "Batch 2",
"3" = "Batch 3", "4" = "Batch 4", "5" = "Batch 5")))
# minimum and maximum cell aspect ratio in the control
AR_min <- morphology_abundance_data %>% ungroup() %>% filter(batch == 1 & day > 0) %>%
select(mean_AR_day) %>% min()
AR_max <- morphology_abundance_data %>% ungroup() %>% filter(batch == 1 & day > 0) %>%
select(mean_AR_day) %>% max()
# plot cell aspect ratio
plot_cell_AR <- ggplot(data = morphology_abundance_data, aes(x = batch_day, y = mean_AR_day,
group = population, col = name)) +
geom_point(size = 1) +
geom_hline(aes(yintercept = AR_min), linetype = "dashed") +
geom_hline(aes(yintercept = AR_max), linetype = "dashed") +
geom_errorbar(aes(ymin = mean_AR_day - se_AR_day, ymax = mean_AR_day + se_AR_day, width = 0.9), size = 0.2) +
geom_line() +
xlab("") +
ylab("Cell shape") +
scale_color_manual(values = c("#d2e0ee", "#efd0d4", "#90b2d5", "#d88b95", "#4d84bc","#c14655", "#2166ac","#b2182b"),
name = "Population and\n temperature",
labels = as_labeller(c("1_20" = "Pop. 1 20°C", "1_38" = "Pop. 1 38°C",
"2_20" = "Pop. 2 20°C", "2_38" = "Pop. 2 38°C",
"3_20" = "Pop. 3 20°C", "3_38" = "Pop. 3 38°C",
"4_20" = "Pop. 4 20°C", "4_38" = "Pop. 4 38°C"))) +
scale_x_continuous(breaks = c(5, 10)) +
theme(text = element_text(size = large_font), legend.text = element_text(size = small_font),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
legend.position = "bottom", strip.text = element_blank(),
plot.margin = unit(c(0, 1.5, 0 ,1.5 ),"cm")) +
facet_wrap(~ batch, ncol = 5, labeller = as_labeller(c("1" = "Batch 1", "2" = "Batch 2",
"3" = "Batch 3", "4" = "Batch 4", "5" = "Batch 5")))
# all plots in one figure
all_time_series <- ggarrange(plot_abundances, plot_cell_area, plot_var_cell_area, plot_cv_cell_area, plot_cell_AR,
nrow = 5, common.legend = T, legend = "bottom", align = "hv",
labels = "AUTO")
# print the plot
all_time_series
# load the package growthrates
library(growthrates)
# select only relevant data for growth model (population, population density and temperature)
# exclude populations 32, 141 and 341, since model will be run separetely for them
data_subset1 <- morphology_abundance_data %>%
group_by(population, lineage, batch, replicate, temperature) %>%
select(population, mean_density, batch_day) %>%
filter(!population %in% c("32", "141", "341"))
# set parameters to start the model (p), lower boundary and upper boundary
# y0 = initial population density, mumax = maximum growth rate, K = carrying capacity, lambda = lag phase
p <- c(y0 = 100, mumax = 1, K = 1000000, lambda = 5)
lower <- c(y0 = 0, mumax = 0, K = 100000, lambda = 0)
upper <- c(y0 = 5000, mumax = 15, K = 5000000, lambda = 8)
# run gompertz model with lag phase estimation for each population separately
subset1_gompertz3 <- all_growthmodels(
mean_density ~ grow_gompertz3(batch_day, parms) | population,
data = data_subset1,
p = p, lower = lower, upper = upper,
log = "y", ncores = 4)
# check the rsquared of the models
# r squared of the model
rsquared(subset1_gompertz3) > 0.9
## 11 12 21 22 31 41 42 131 132 142 151 152
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 231 232 241 242 251 252 331 332 342 351 352 431
## TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 432 441 442 451 452
## TRUE TRUE TRUE TRUE TRUE
rsquared(subset1_gompertz3) # all above 0.90, expect pop. 241, with 0.80
## 11 12 21 22 31 41 42
## 0.9883432 0.9178892 0.9933680 0.9998310 0.9937436 0.9809305 0.9981437
## 131 132 142 151 152 231 232
## 0.9998194 0.9998437 0.9998797 0.9914330 0.9998212 0.9999818 0.9999183
## 241 242 251 252 331 332 342
## 0.8075650 0.9992203 0.9998438 0.9998244 0.9998084 0.9539876 0.9999996
## 351 352 431 432 441 442 451
## 0.9477612 0.9993824 0.9992420 0.9996416 0.9997832 0.9999890 0.9997705
## 452
## 0.9994084
# coefficients of the model
coef(subset1_gompertz3)
## y0 mumax K lambda
## 11 4999.99925 1.197897 1129221.8 5.426150e-06
## 12 490.11662 3.939040 156097.9 7.999999e+00
## 21 2822.39743 1.342365 1261845.6 1.789159e-06
## 22 762.82078 1.764727 584705.5 7.973149e+00
## 31 627.35002 1.906846 1173617.9 2.705005e-06
## 41 4999.99917 1.173631 1373269.9 4.339097e-06
## 42 1736.92242 1.354624 100000.0 7.647302e+00
## 131 1654.64387 1.676791 1081410.9 7.102917e+00
## 132 1589.40401 1.251949 3011394.3 6.268673e+00
## 142 1452.04283 1.144168 4999513.5 4.818932e-01
## 151 19.72001 4.887823 248032.6 4.603856e-01
## 152 1335.03526 1.412242 2279297.1 2.774922e-01
## 231 1426.42616 2.897612 504515.5 7.408923e+00
## 232 924.84138 2.550602 4997991.7 7.924815e+00
## 241 3241.76987 14.999383 226100.6 2.869473e+00
## 242 4999.99967 2.225997 743961.4 2.876162e+00
## 251 399.02617 2.094700 4685750.6 1.448947e+00
## 252 855.05594 2.145062 924045.8 9.043278e-01
## 331 2093.40965 1.539841 4456116.9 3.773430e+00
## 332 4313.96423 14.999743 264803.7 4.813319e+00
## 342 806.37213 1.727175 1957003.6 6.825973e-01
## 351 2024.29552 14.760898 201141.4 1.765625e+00
## 352 1050.07698 1.509381 1850540.4 2.545231e-01
## 431 107.09779 2.026098 274606.3 3.168215e+00
## 432 1341.19534 2.204486 137174.9 4.250084e+00
## 441 1618.63891 1.823084 167147.9 2.417790e+00
## 442 1558.57363 2.097457 767992.8 1.765452e+00
## 451 2814.32870 3.254367 264303.9 1.972106e+00
## 452 583.12852 2.192975 808572.3 7.190942e-01
# plot the model, y in log scale
par(mfrow = c(8, 4))
par(mar = c(2.5, 4, 2, 1))
plot(subset1_gompertz3)
plot(subset1_gompertz3, log = "y")
# populations 32, 141, 341 presented suboptimal fittings with the above parameters
# model these populations separately
# population 32, remove last day since population is already crashing
pop32 <- morphology_abundance_data %>%
group_by(population, lineage, batch, replicate, temperature) %>%
select(population, mean_density, batch_day) %>%
filter(population == 32 & batch_day < 12)
## Adding missing grouping variables: `lineage`, `batch`, `replicate`, `temperature`
# same parameters as previous populations
pop32_gompertz <- fit_growthmodel(FUN = grow_gompertz3, p = p,
pop32$batch_day, pop32$mean_density,
lower = lower, upper = upper)
# # population 141
pop141 <- morphology_abundance_data %>%
group_by(population, lineage, batch, replicate, temperature) %>%
select(population, mean_density, batch_day) %>%
filter(population == 141 & batch_day < 6)
## Adding missing grouping variables: `lineage`, `batch`, `replicate`, `temperature`
# use new parameters
p2 <- c(y0 = 500, mumax = 1, K = 1000000, lambda = 0.1)
lower2 <- c(y0 = 0, mumax = 0, K = 100000, lambda = 0)
upper2 <- c(y0 = 10000, mumax = 30, K = 5000000, lambda = 4)
pop141_gompertz <- fit_growthmodel(FUN = grow_gompertz3, p = p,
pop141$batch_day, pop141$mean_density,
lower = lower, upper = upper)
# population 341, remove final day that population is crashing
pop341 <- morphology_abundance_data %>%
group_by(population, lineage, batch, replicate, temperature) %>%
select(population, mean_density, batch_day) %>%
filter(population == 341 & batch_day <= 5)
## Adding missing grouping variables: `lineage`, `batch`, `replicate`, `temperature`
# rsquared of the model is very low, manually estimate growth rate and lag phase
pop341_day0 <- pop341 %>% filter(batch_day == 0)
pop341_day3 <- pop341 %>% filter(batch_day == 3)
# calculate maximum growth rate as the log difference in population density divided by the time difference
mumax <- (log10(pop341_day3$mean_density) - log10(pop341_day0$mean_density))/(pop341_day3$batch_day - pop341_day0$batch_day)
lambda <- 0 # define lag phase as 0
y0 <- NA
K <- NA
population <- 341
pop341_manual_growth <- tibble(mumax, lambda, y0, K, population)
remove(mumax, lambda, y0, K, population)
# pop 32
rsquared(pop32_gompertz)
## [1] 0.8966222
coef(pop32_gompertz)
## y0 mumax K lambda
## 2.969880e+03 7.143717e+00 1.681196e+05 7.999997e+00
par(mfrow = c(1,2))
plot(pop32_gompertz)
plot(pop32_gompertz, log = "y")
# pop 141
rsquared(pop141_gompertz)
## [1] 0.9996588
coef(pop141_gompertz)
## y0 mumax K lambda
## 2.047749e+03 3.552567e+00 2.361080e+05 2.279727e+00
par(mfrow = c(1,2))
plot(pop141_gompertz)
plot(pop141_gompertz, log = "y")
# try to get the model predictions for all gompertz models
# new time variable with the maximum batch_day across all populations (12)
time <- data.frame(time = seq(from = 0, to = 12, by = 1))
# predict function for each problematic population, add name of the population to the tibble
pop32_growth_prediction <- pop32_gompertz %>% predict(newdata = time) %>% as_tibble() %>% dplyr::rename(batch_day = time, mean_density_predict = y) %>% mutate(population = "32")
pop141_growth_prediction <- pop141_gompertz %>% predict(newdata = time) %>% as_tibble() %>% dplyr::rename(batch_day = time, mean_density_predict = y) %>% mutate(population = "141")
# predict function for each the remaining populations
subset1_growth_prediction <- predict(subset1_gompertz3, newdata = time)
names(subset1_growth_prediction)
## [1] "11" "12" "21" "22" "31" "41" "42" "131" "132" "142" "151"
## [12] "152" "231" "232" "241" "242" "251" "252" "331" "332" "342" "351"
## [23] "352" "431" "432" "441" "442" "451" "452"
# transform list into data frame with population as a column
subset1_growth_prediction <- map_df(subset1_growth_prediction, ~as.data.frame(.x), .id = "population")
# rename the columns
subset1_growth_prediction <- subset1_growth_prediction %>% dplyr::rename(batch_day = time, mean_density_predict = y)
# merge all models
prediction_gompertz <- bind_rows(subset1_growth_prediction, pop32_growth_prediction, pop141_growth_prediction)
prediction_gompertz$population <- as.numeric(prediction_gompertz$population)
morphology_abundance_data_gompertz <- morphology_abundance_data %>% left_join(prediction_gompertz)
# create a dataset with the coeficients of the model
coef_gompertz <- coef(subset1_gompertz3) %>% as_tibble()
coef_gompertz <- cbind(coef_gompertz, population = as.double(names(subset1_gompertz3)))
coef_gompertz <- coef_gompertz %>% rbind(c(coef(pop32_gompertz), 32))
coef_gompertz <- coef_gompertz %>% rbind(c(coef(pop141_gompertz), 141))
coef_gompertz <- coef_gompertz %>% rbind(pop341_manual_growth)
# select metadata from the original dataset
metadata <- morphology_abundance_data %>% select(population, lineage, batch, replicate, temperature) %>% distinct()
# add metadata to the growth model coeficients data
coef_gompertz <- coef_gompertz %>% left_join(metadata)
## Joining, by = "population"
# add name column for plotting (lineage + temperature)
coef_gompertz <- coef_gompertz %>% mutate(name = paste0(lineage, "_", temperature))
# define font sizes
small_font <- 10
large_font <- 13
# plot the lag phase length (lambda)
plot_lambda <- ggplot(coef_gompertz, aes(x = batch, y = lambda, fill = name, group = population)) +
geom_col(position = "dodge") +
xlab("") +
ylab("Lag phase (days)") +
scale_fill_manual(values = c("#d2e0ee", "#efd0d4", "#90b2d5", "#d88b95", "#4d84bc","#c14655", "#2166ac","#b2182b"),
name = "Population and\n temperature",
labels = as_labeller(c("1_20" = "Pop. 1 20°C", "1_38" = "Pop. 1 38°C",
"2_20" = "Pop. 2 20°C", "2_38" = "Pop. 2 38°C",
"3_20" = "Pop. 3 20°C", "3_38" = "Pop. 3 38°C",
"4_20" = "Pop. 4 20°C", "4_38" = "Pop. 4 38°C"))) +
guides(fill = guide_legend(title.theme = element_text(size = small_font), label.theme = element_text(size = small_font))) +
theme(text = element_text(size = large_font),
panel.border = element_rect(colour = "black", fill = NA),
panel.background = element_blank(),
legend.position = "right",
#plot.margin = margin(t = 10, r = 10, b = 10, l = 10),
axis.title.y = element_text(margin = margin(t = 0, r = 10, b = 0, l = 0)),
axis.text.x = element_blank(), axis.ticks.x = element_blank()) + #,
facet_grid(~batch, scales = "free_x",
labeller = as_labeller(c("1" = "Batch 1", "2" = "Batch 2",
"3" = "Batch 3", "4" = "Batch 4", "5" = "Batch 5")))
# plot the maximum growth rate (mumax)
plot_mumax <- ggplot(coef_gompertz, aes(x = batch, y = mumax, fill = name, group = population)) +
geom_col(position = "dodge") +
xlab("") +
ylab(bquote("Max. growth rate"~(day^-1))) +
scale_fill_manual(values = c("#d2e0ee", "#efd0d4", "#90b2d5", "#d88b95", "#4d84bc","#c14655", "#2166ac","#b2182b"),
name = "Population and\n temperature", guide = F,
labels = as_labeller(c("1_20" = "Pop. 1 20°C", "1_38" = "Pop. 1 38°C",
"2_20" = "Pop. 2 20°C", "2_38" = "Pop. 2 38°C",
"3_20" = "Pop. 3 20°C", "3_38" = "Pop. 3 38°C",
"4_20" = "Pop. 4 20°C", "4_38" = "Pop. 4 38°C"))) +
#guides(fill = guide_legend(title.theme = element_text(size = small_font), label.theme = element_text(size = small_font))) +
theme(text = element_text(size = large_font),
panel.border = element_rect(colour = "black", fill = NA),
panel.background = element_blank(),
legend.position = "bottom",
#plot.margin = margin(t = 10, r = 10, b = 10, l = 10),
axis.title.y = element_text(margin = margin(t = 0, r = 10, b = 0, l = 0)),
axis.text.x = element_blank(), axis.ticks.x = element_blank()) + #,
facet_grid(~batch, scales = "free_x",
labeller = as_labeller(c("1" = "Batch 1", "2" = "Batch 2",
"3" = "Batch 3", "4" = "Batch 4", "5" = "Batch 5")))
gompertz_fit_plot <- ggplot(data = morphology_abundance_data_gompertz) +
geom_point(aes(x = batch_day, y = mean_density, group = population, colour = name), size = 1.5) +
geom_line(aes(x = batch_day, y = mean_density_predict, group = population, colour = name)) +
scale_y_log10() +
ylab("Population density (cells/ml)\n") +
xlab("Day\n") +
scale_color_manual(values = c("#d2e0ee", "#efd0d4", "#90b2d5", "#d88b95", "#4d84bc","#c14655", "#2166ac","#b2182b", "red", "blue"),
name = "Population and\n temperature", guide = F,
labels = as_labeller(c("1_20" = "Pop. 1 20°C", "1_38" = "Pop. 1 38°C",
"2_20" = "Pop. 2 20°C", "2_38" = "Pop. 2 38°C",
"3_20" = "Pop. 3 20°C", "3_38" = "Pop. 3 38°C",
"4_20" = "Pop. 4 20°C", "4_38" = "Pop. 4 38°C"))) +
scale_x_continuous(breaks = c(5, 10)) +
theme(text = element_text(size = large_font),
panel.background = element_blank(), panel.border = element_rect(colour = "black", fill = NA),
axis.line = element_line(colour = "black"),
legend.position = "bottom",
plot.margin = unit(c(0, 1.5, 0 , 1.5 ),"cm")) +
facet_grid(paste0("Batch ", batch) ~ paste0("Population ", lineage))
# layout of the final figure
layout <- "
AA#
AA#
AA#
BBD
CCD"
# create final figure with three facets
plot_growth <- gompertz_fit_plot + plot_mumax + plot_lambda + guide_area() + plot_layout(guides = "collect", design = layout) + plot_annotation(tag_levels = 'A')
# print figure
plot_growth
# create variable merging batch, temperature and replicate; select only used variables
growth_response_gompertz <- coef_gompertz %>% ungroup() %>%
mutate(batch_temp = paste0("b", batch, "_", temperature, "_", replicate)) %>%
select(-population, -batch, - replicate, - temperature, -name)
# calculate % change in comparison to the same population lineage in batch 1
growth_response_gompertz <- growth_response_gompertz %>% select(-y0, -K) %>% pivot_longer(cols = 1:2, names_to = "data_type")
growth_response_gompertz <- growth_response_gompertz %>% group_by(lineage, data_type) %>% pivot_wider(names_from = batch_temp, values_from = value)
# calculate the difference in maximum growth rate and in lag phase duration in absolute values, not the % change since numbers are too large and interpretation is difficult
growth_response_gompertz <- growth_response_gompertz %>%
mutate(week2 = b2_38_NA - b1_20_NA,
week3_1 = b3_38_1 - b1_20_NA,
week3_2 = b3_38_2 - b1_20_NA,
week4 = b4_38_1 - b1_20_NA,
week5 = b5_38_1 - b1_20_NA,
week4_cold = b4_20_2 - b1_20_NA,
week5_cold = b5_20_2 - b1_20_NA)
# tidy the data
growth_response_gompertz <- growth_response_gompertz %>% select(1, 2, 11:17) %>%
gather(key = "comparison", value = "percent_batch1", 3:9)
# repeat the same calculations with the morphological traits
# calculate mean cell area and cell aspect ratio in each batch
morphology_data_batch <- morphology_abundance_data %>% drop_na(mean_area_day) %>%
group_by(population, lineage, batch, replicate, temperature) %>%
summarise(mean_area_batch = mean(mean_area_day), cv_area_batch = mean(cv_area_day), var_area_batch = mean(var_area_day),
mean_AR_batch = mean(mean_AR_day))
# create variable combining batch, temperature and replicate; select only used variables
morphology_response <- morphology_data_batch %>% ungroup() %>%
mutate(batch_temp = paste0("b",batch, "_", temperature, "_", replicate)) %>%
select(lineage, batch_temp, mean_area_batch, mean_AR_batch, cv_area_batch, var_area_batch)
morphology_response <- morphology_response %>% gather(key = "data_type", value = "value", 3:6)
morphology_response <- morphology_response %>% group_by(lineage, data_type) %>% spread(key = batch_temp, value = value)
# calculate % change in comparison to the same population lineage in batch 1
morphology_response <- morphology_response %>% group_by(lineage, data_type) %>%
mutate(week2 = ((b2_38_NA - b1_20_NA) * 100) / b1_20_NA,
week3_1 = ((b3_38_1 - b1_20_NA) * 100) / b1_20_NA,
week3_2 = ((b3_38_2 - b1_20_NA) * 100) / b1_20_NA,
week4 = ((b4_38_1 - b1_20_NA) * 100) / b1_20_NA,
week5 = ((b5_38_1 - b1_20_NA) * 100) / b1_20_NA,
week4_cold = ((b4_20_2 - b1_20_NA) * 100) / b1_20_NA,
week5_cold = ((b5_20_2 - b1_20_NA) * 100) / b1_20_NA)
morphology_response <- morphology_response %>% select(1, 2, 11:17) %>%
gather(key = "comparison", value = "percent_batch1", 3:9)
# merge growth and lag phase with the morphology data
response_data <- rbind.data.frame(growth_response_gompertz, morphology_response)
# add batch information
response_data <- response_data %>% ungroup() %>% mutate(batch = case_when(
comparison == "week2" ~ 2,
comparison == "week3_1" ~ 3,
comparison == "week3_2" ~ 3,
comparison == "week4" ~ 4,
comparison == "week5" ~ 5,
comparison == "week4_cold" ~ 4,
comparison == "week5_cold" ~ 5
))
# add temperature information
response_data <- response_data %>% mutate(temperature = case_when(
comparison == "week2" ~ 38,
comparison == "week3_1" ~ 38,
comparison == "week3_2" ~ 38,
comparison == "week4" ~ 38,
comparison == "week5" ~ 38,
comparison == "week4_cold" ~ 20,
comparison == "week5_cold" ~ 20
))
# add variable for lineage and temperature color
response_data <- response_data %>% mutate(name = paste0(lineage, "_", temperature))
# change data type into factor and change levels for a nicer looking plot
response_data$data_type <- factor(response_data$data_type, levels = c("mumax", "lambda",
"mean_area_batch", "cv_area_batch", "var_area_batch",
"mean_AR_batch"))
# transform lineage into categorical variable
response_data$lineage <- as.character(response_data$lineage)
# look at the data
ggplot(data = response_data, aes(x = batch, y = percent_batch1, color = name)) +
geom_point(alpha = 0.8, size = 2) +
geom_hline(yintercept = 0, linetype = "dashed") +
scale_color_manual(values = c("#d2e0ee", "#efd0d4", "#90b2d5", "#d88b95", "#4d84bc","#c14655", "#2166ac","#b2182b", "#878787"),
name = "Population and\n temperature",
labels = as_labeller(c("1_20" = "Pop. 1 20°C", "1_38" = "Pop. 1 38°C",
"2_20" = "Pop. 2 20°C", "2_38" = "Pop. 2 38°C",
"3_20" = "Pop. 3 20°C", "3_38" = "Pop. 3 38°C",
"4_20" = "Pop. 4 20°C", "4_38" = "Pop. 4 38°C"))) +
ylab("Change compared to control (%)") +
xlab("Batch culture") +
guides(fill = guide_legend(title.theme = element_text(size = 15),
label.theme = element_text(size = 15))) +
theme(text = element_text(size = 15),
panel.border = element_rect(colour = "black", fill = NA),
panel.background = element_blank(),
legend.position = "bottom") +
facet_wrap(~data_type, ncol = 3, labeller = as_labeller(c("mumax" = "Growth rate",
"lambda" = "Lag phase",
"cv_area_batch" = "CV cell size",
"mean_AR_batch" = "Cell shape",
"mean_area_batch" = "Cell size",
"var_area_batch" = "Variance cell size")), scales = "free_y")
# load the libraries
library(MCMCglmm)
library(coda)
library(purrr)
set.seed(1)
# select the populations that experienced 38 °C
response_data_38 <- response_data %>% filter(temperature == 38)
# create variable for time difference since warming started
response_data_38 <- response_data_38 %>% mutate(batches_38 = batch - 2)
# create a list with the prior information, using an uninformative prior
prior_1 <- list(R = list(V = 1, nu = 0.002),
G = list(G1 = list(V = 1, nu = 0.002)))
# run one separate model for each of the 5 response variables
models_38 <- response_data_38 %>% group_by(data_type) %>% nest() %>%
mutate(mixed_model = map(data, ~MCMCglmm(fixed = percent_batch1 ~ batches_38,
random = ~lineage,
data = .,
family = "gaussian",
prior = prior_1,
nitt = 2000000,
burnin = 30000,
thin = 1000, verbose = F)))
# create new data
new_data <- data.frame(expand.grid(lineage = c("1"), batches_38 = c(0:3),
percent_batch1 = 0))
# get confidence interval from each model
models_38 <- models_38 %>% mutate(predictions = map(mixed_model, ~predict.MCMCglmm(object = .,
newdata = new_data,
interval = "confidence")))
models_38 <- models_38 %>% mutate(predictions = map(.x = predictions, ~as_tibble(.)))
# merge new data and model fit for plotting
models_38 <- models_38 %>% mutate(predictions = map(.x = predictions, ~cbind(new_data, .)))
# store models in a separate object
model_38_predictions <- unnest(models_38, predictions)
# plot the data and the model
response_38_plot <- ggplot(data = response_data_38, aes(x = batch, y = percent_batch1, color = name)) +
geom_point(aes(color = as.factor(lineage)), alpha = 0.8, size = 2) +
geom_smooth(data = model_38_predictions, aes(x = batches_38 + 2, y = fit, ymin = lwr, ymax = upr), colour = "black", stat = "identity") +
geom_hline(yintercept = 0, linetype = "dashed") +
#scale_y_continuous(limits = c(-95, 95)) +
xlab("Batch") +
ylab("Change in comparison to control (%)") +
scale_color_manual(values = c("#efd0d4", "#d88b95","#c14655","#b2182b"),
name = "Population", labels = as_labeller(c( "1" = "Pop. 1",
"2" = "Pop. 2",
"3" = "Pop. 3",
"4" = "Pop. 4"))) +
theme(text = element_text(size = 15), plot.title = element_text(hjust = 0.5, margin = unit(c(0.5,0.5,0.5,0.5), "lines")),
panel.border = element_rect(colour = "black", fill = NA),
panel.background = element_blank(),
legend.position = "bottom") +
facet_wrap(~data_type, ncol = 3, labeller = as_labeller(c("mumax" = "Growth rate",
"lambda" = "Lag phase",
"cv_area_batch" = "CV cell size",
"mean_AR_batch" = "Cell shape",
"mean_area_batch" = "Cell size",
"var_area_batch" = "Variance cell size")), scales = "free_y")
response_38_plot
# select the populations that experienced 20 °C
response_data_20 <- response_data %>% filter(temperature == 20)
# create variable for time difference since return to 20 °C started
response_data_20 <- response_data_20 %>% mutate(batches_20 = batch - 4)
# create a list with the prior information, using an uninformative prior
# same prior as 38°C models
prior_1_new <- list(R = list(V = 1, nu = 0.002),
G = list(G1 = list(V = 1, nu = 0.002)))
# run one separate model for each of the 5 response variables
models_20 <- response_data_20 %>% group_by(data_type) %>% nest() %>%
mutate(mixed_model = map(data, ~MCMCglmm(fixed = percent_batch1 ~ batches_20,
random = ~lineage,
data = .,
family = "gaussian",
prior = prior_1_new,
nitt = 2000000,
burnin = 30000,
thin = 1000,
verbose = F)))
# create new data
new_data_20 <- data.frame(expand.grid(lineage = c("1"), batches_20 = c(0, 1),
percent_batch1 = 0))
# get confidence interval from each model
models_20 <- models_20 %>% mutate(predictions = map(mixed_model, ~predict.MCMCglmm(object = .,
newdata = new_data_20,
interval = "confidence")))
models_20 <- models_20 %>% mutate(predictions = map(.x = predictions, ~as_tibble(.)))
# merge new data and model fit for plotting
models_20 <- models_20 %>% mutate(predictions = map(.x = predictions, ~cbind(new_data_20, .)))
# store models in a separate object
model_20_predictions <- unnest(models_20, predictions)
# plot the data and the model
response_20_plot <- ggplot(data = response_data_20, aes(x = batch, y = percent_batch1, color = name)) +
geom_point(aes(color = as.factor(lineage)), alpha = 0.8, size = 2) +
geom_smooth(data = model_20_predictions,
aes(x = batches_20 + 4, y = fit, ymin = lwr, ymax = upr),
colour = "black", stat = "identity") +
geom_hline(yintercept = 0, linetype = "dashed") +
scale_x_continuous(breaks = c(4, 5)) +
xlab("Batch") +
ylab("Change in comparison to control (%)") +
scale_color_manual(values = c("#d2e0ee", "#90b2d5", "#4d84bc", "#2166ac"),
name = "Population",
labels = as_labeller(c( "1" = "Pop. 1",
"2" = "Pop. 2",
"3" = "Pop. 3",
"4" = "Pop. 4"))) +
theme(text = element_text(size = 15), plot.title = element_text(hjust = 0.5, margin = unit(c(0.5,0.5,0.5,0.5), "lines")),
panel.border = element_rect(colour = "black", fill = NA),
panel.background = element_blank(),
legend.position = "bottom") +
facet_wrap(~data_type, ncol = 3, labeller = as_labeller(c("mumax" = "Growth rate",
"lambda" = "Lag phase",
"cv_area_batch" = "CV cell size",
"mean_AR_batch" = "Cell shape",
"mean_area_batch" = "Cell size",
"var_area_batch" = "Variance cell size")), scales = "free_y")
response_20_plot
# plot both models in the same plot with all the response data
response_plot <- ggplot(data = filter(response_data, data_type != "K"), aes(x = batch, y = percent_batch1, color = name)) +
geom_point(alpha = 0.8, size = 2) +
geom_hline(yintercept = 0, linetype = "dashed") +
scale_color_manual(values = c("#d2e0ee", "#efd0d4", "#90b2d5", "#d88b95", "#4d84bc","#c14655", "#2166ac","#b2182b"),
name = "Population and\n temperature",
labels = as_labeller(c("1_20" = "Pop. 1 20°C", "1_38" = "Pop. 1 38°C",
"2_20" = "Pop. 2 20°C", "2_38" = "Pop. 2 38°C",
"3_20" = "Pop. 3 20°C", "3_38" = "Pop. 3 38°C",
"4_20" = "Pop. 4 20°C", "4_38" = "Pop. 4 38°C"))) +
geom_smooth(data = model_20_predictions %>% filter(data_type != "K"),
aes(x = batches_20 + 4, y = fit, ymin = lwr, ymax = upr), alpha = 0.2, colour = "#82a8d0", stat = "identity") +
geom_smooth(data = model_38_predictions %>% filter(data_type != "K"),
aes(x = batches_38 + 2, y = fit, ymin = lwr, ymax = upr), alpha = 0.2, colour = "#d37d88", stat = "identity") +
geom_hline(yintercept = 0, linetype = "dashed") +
xlab("Batch") +
ylab("Change in comparison to control (%)") +
theme(text = element_text(size = 15),
panel.border = element_rect(colour = "black", fill = NA),
panel.background = element_blank(),
legend.position = "bottom",
plot.margin = unit(c(1.5, 1.5, 1.5 ,1.5 ),"cm")) +
facet_wrap(~data_type, ncol = 3, labeller = as_labeller(c("mumax" = "Growth rate",
"lambda" = "Lag phase",
"cv_area_batch" = "CV cell size",
"mean_AR_batch" = "Cell shape",
"mean_area_batch" = "Cell size",
"var_area_batch" = "Variance cell size")), scales = "free_y")
response_plot
model_38_predictions %>% group_by(data_type) %>% select(-percent_batch1, -lineage) %>% filter(batches_38 == 3)
## # A tibble: 6 x 7
## # Groups: data_type [6]
## data_type data mixed_model batches_38 fit lwr upr
## <fct> <list<df[,7> <list> <int> <dbl> <dbl> <dbl>
## 1 mumax [20 × 7] <MCMCglmm> 3 4.72 0.468 9.08
## 2 lambda [20 × 7] <MCMCglmm> 3 0.756 -0.614 2.15
## 3 cv_area_batch [20 × 7] <MCMCglmm> 3 32.9 5.75 56.6
## 4 mean_AR_batch [20 × 7] <MCMCglmm> 3 -16.5 -19.0 -14.1
## 5 mean_area_bat… [20 × 7] <MCMCglmm> 3 -30.2 -49.8 -11.2
## 6 var_area_batch [20 × 7] <MCMCglmm> 3 -17.4 -39.5 4.79
model_38_predictions %>% group_by(data_type) %>% select(-percent_batch1, -lineage) %>% filter(batches_38 == 0)
## # A tibble: 6 x 7
## # Groups: data_type [6]
## data_type data mixed_model batches_38 fit lwr upr
## <fct> <list<df[,7]>> <list> <int> <dbl> <dbl> <dbl>
## 1 mumax [20 × 7] <MCMCglmm> 0 1.58 -2.30 5.68
## 2 lambda [20 × 7] <MCMCglmm> 0 7.72 6.46 9.09
## 3 cv_area_batch [20 × 7] <MCMCglmm> 0 43.5 19.6 69.2
## 4 mean_AR_batch [20 × 7] <MCMCglmm> 0 -10.7 -12.9 -8.38
## 5 mean_area_batch [20 × 7] <MCMCglmm> 0 -27.3 -47.5 -9.27
## 6 var_area_batch [20 × 7] <MCMCglmm> 0 13.5 -6.14 36.0
model_20_predictions %>% group_by(data_type) %>% select(-percent_batch1, -lineage) %>% filter(batches_20 == 1)
## # A tibble: 6 x 7
## # Groups: data_type [6]
## data_type data mixed_model batches_20 fit lwr upr
## <fct> <list<df[,> <list> <dbl> <dbl> <dbl> <dbl>
## 1 mumax [8 × 7] <MCMCglmm> 1 0.442 -0.438 1.47
## 2 lambda [8 × 7] <MCMCglmm> 1 0.557 -0.607 1.64
## 3 cv_area_batch [8 × 7] <MCMCglmm> 1 -10.3 -13.6 -6.62
## 4 mean_AR_batch [8 × 7] <MCMCglmm> 1 -5.16 -10.6 0.0956
## 5 mean_area_ba… [8 × 7] <MCMCglmm> 1 -11.1 -16.6 -5.01
## 6 var_area_bat… [8 × 7] <MCMCglmm> 1 -36.9 -43.8 -30.5